---
title: "Constitutions"
author: "Sally Sharif"
date: "2024-02-10"
output: html_document
---

```{r}
library(dplyr)
library(haven)
library(table1)
library(corrtable)
library(car)
library(stargazer)
library(rio)
library(mice)
library(mitools)

## Load the data
data <- read_dta("rebel_constitutions.dta")

# Rename variables for better display
data <- data %>% rename(
  country_code = cowcode,
  year = year,
  private_property = proprght,
  expropriation_gpp = expcond_7,
  expropriation_pay = explim_2,
  expropriation_legal = explim_4,
  rebel_regime = rebel_regime,
  population = WB_population,
  gdp_per_capita = WB_gdp_pc,
  gdp_growth = WB_gdp_growth,
  oil_percent = WB_oil_percent,
  ethnic_fractionalization = ethnic,
  religious_fractionalization = religious,
  coups = coups,
  democracy = democracy
)

## Descriptive statistics
data_desc <- data %>%
  mutate(
    expropriation_pay = as.factor(expropriation_pay),
    expropriation_gpp = as.factor(expropriation_gpp),
    expropriation_legal = as.factor(expropriation_legal),
    private_property = as.factor(private_property),
    rebel_regime = as.factor(rebel_regime),
    democracy = as.factor(democracy)
  )

# Appendix Tables 2-3
table1(~ rebel_regime + private_property + expropriation_gpp + expropriation_pay + expropriation_legal + population + gdp_per_capita + gdp_growth + oil_percent + ethnic_fractionalization + religious_fractionalization + coups + democracy, data = data_desc)

## Correlations
data$private_property[data$private_property %in% c(90, 96)] <- NA
data$expropriation_gpp[data$expropriation_gpp == 99] <- NA
data$expropriation_pay[data$expropriation_pay == 99] <- NA
data$expropriation_legal[data$expropriation_legal == 99] <- NA

# Appendix Table 4
correlation_matrix(data[, c("rebel_regime", "private_property", "expropriation_gpp", "expropriation_pay", "expropriation_legal", "population", "gdp_per_capita", "gdp_growth", "oil_percent", "ethnic_fractionalization", "religious_fractionalization", "coups", "democracy")], digits = 2, use="lower", replace_diagonal = TRUE)
```

```{r}
# Logistic regression models
# Manuscript Tables 2-3

run_logistic_models <- function(dv) {
  data_clean <- data %>% filter(get(dv) %in% c(0, 1))
  
  model1 <- glm(as.formula(paste(dv, "~ rebel_regime")), data = data_clean, family = binomial)
  model2 <- glm(as.formula(paste(dv, "~ rebel_regime + log(population) + log(gdp_per_capita) + gdp_growth + oil_percent + ethnic_fractionalization + religious_fractionalization + coups")), data = data_clean, family = binomial)
  model3 <- glm(as.formula(paste(dv, "~ rebel_regime + log(population) + log(gdp_per_capita) + gdp_growth + oil_percent + ethnic_fractionalization + religious_fractionalization + coups + factor(year)")), data = data_clean, family = binomial)
  
  list(
    summary1 = summary(model1),
    summary2 = summary(model2),
    summary3 = summary(model3),
    vif2 = vif(model2)
  )
}

# Run models for each dependent variable
results_private_property_logit <- run_logistic_models("private_property")
results_expropriation_gpp_logit <- run_logistic_models("expropriation_gpp")
results_expropriation_pay_logit <- run_logistic_models("expropriation_pay")
results_expropriation_legal_logit <- run_logistic_models("expropriation_legal")

# Print results separately
print(results_private_property_logit)
print(results_expropriation_gpp_logit)
print(results_expropriation_pay_logit)
print(results_expropriation_legal_logit)
```

```{r}
# Logistic regression models with alternative controls (democracy instead of oil)
# Appendix Tables 6-7

run_models_alt_controls <- function(dv) {
  data_clean <- data %>% filter(get(dv) %in% c(0, 1))
  
  model1 <- glm(as.formula(paste(dv, "~ rebel_regime")), data = data_clean, family = binomial)
  model2 <- glm(as.formula(paste(dv, "~ rebel_regime + log(population) + log(gdp_per_capita) + gdp_growth + democracy + ethnic_fractionalization + religious_fractionalization + coups")), data = data_clean, family = binomial)
  model3 <- glm(as.formula(paste(dv, "~ rebel_regime + log(population) + log(gdp_per_capita) + gdp_growth + democracy + ethnic_fractionalization + religious_fractionalization + coups + factor(year)")), data = data_clean, family = binomial)
  
  list(
    summary1 = summary(model1),
    summary2 = summary(model2),
    summary3 = summary(model3),
    vif2 = vif(model2)
  )
}

# Run models for each dependent variable
results_private_property_alt <- run_models_alt_controls("private_property")
results_expropriation_gpp_alt <- run_models_alt_controls("expropriation_gpp")
results_expropriation_pay_alt <- run_models_alt_controls("expropriation_pay")
results_expropriation_legal_alt <- run_models_alt_controls("expropriation_legal")

# Print results separately
print(results_private_property_alt)
print(results_expropriation_gpp_alt)
print(results_expropriation_pay_alt)
print(results_expropriation_legal_alt)
```

```{r}
## Logistic regression models for countries with at least one year of rebel victory
data_subset <- data %>%
  group_by(country_code) %>%
  filter(any(rebel_regime == 1)) %>%
  ungroup()

# Appendix Tables 8-9

run_models <- function(dv) {
  data_clean <- data_subset %>% filter(get(dv) %in% c(0, 1))
  
  model1 <- glm(as.formula(paste(dv, "~ rebel_regime")), data = data_clean, family = binomial)
  model2 <- glm(as.formula(paste(dv, "~ rebel_regime + log(population) + log(gdp_per_capita) + gdp_growth + oil_percent + ethnic_fractionalization + religious_fractionalization + coups")), data = data_clean, family = binomial)
  model3 <- glm(as.formula(paste(dv, "~ rebel_regime + log(population) + log(gdp_per_capita) + gdp_growth + oil_percent + ethnic_fractionalization + religious_fractionalization + coups + factor(year)")), data = data_clean, family = binomial)
  
  list(
    summary1 = summary(model1),
    summary2 = summary(model2),
    summary3 = summary(model3),
    vif2 = vif(model2)
  )
}

# Run models for each dependent variable
results_private_property <- run_models("private_property")
results_expropriation_gpp <- run_models("expropriation_gpp")
results_expropriation_pay <- run_models("expropriation_pay")
results_expropriation_legal <- run_models("expropriation_legal")

# Print results separately
print(results_private_property)
print(results_expropriation_gpp)
print(results_expropriation_pay)
print(results_expropriation_legal)
```

```{r}
## Multiple imputations
# Select relevant variables for imputation
selected_data <- data %>% select(rebel_regime, private_property, expropriation_gpp, expropriation_pay, expropriation_legal, population, gdp_per_capita, gdp_growth, oil_percent, ethnic_fractionalization, religious_fractionalization, coups)

# Show missing pattern
# Appendix Figure 1
selected_data2 <- data %>% select(rebel_regime, private_property, expropriation_gpp, expropriation_pay, expropriation_legal)
md.pattern(selected_data2, rotate.names = TRUE)

# Appendix Table 10
round(flux(selected_data),3)

# Impute missing values
pt.mice <- mice(selected_data, printFlag = FALSE, m = 100, maxit = 0, method = "cart")
densityplot(pt.mice, layout = c(4,4))

# Run multiple imputation models
pt.mice.dat <- list()
pt.mice.dat$imputations <- lapply(1:pt.mice$m, function(x) complete(pt.mice, x))

# Run multiple imputation models
mice_models <- function(var) {
  lapply(1:pt.mice$m, function(x)
    glm(as.formula(paste(var, "~ rebel_regime + log(population) + log(gdp_per_capita) + gdp_growth + oil_percent + 
                         ethnic_fractionalization + religious_fractionalization + coups")),
        data = complete(pt.mice, x), family = binomial))
}

mice_pooled <- function(models) MIcombine(models) %>% summary()

# Appendix Table 11
mice_pooled(mice_models("private_property"))
mice_pooled(mice_models("expropriation_gpp"))
mice_pooled(mice_models("expropriation_pay"))
mice_pooled(mice_models("expropriation_legal"))
```
